library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsdados <- read_rds("../data/geomorfologia.rds")
glimpse(dados)## Rows: 106
## Columns: 22
## $ sup <chr> "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I…
## $ solo <chr> "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "L…
## $ amostra <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
## $ x <dbl> 0, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 32…
## $ amg <dbl> 0.2, 0.1, 0.7, 0.4, 0.4, 0.4, 1.2, 0.8, 1.1, 1.2, 0.1, 0.2, 0.…
## $ ag <dbl> 3.72, 4.27, 5.00, 3.80, 3.10, 3.80, 3.60, 4.70, 4.50, 5.90, 5.…
## $ am <dbl> 20.4, 22.6, 22.7, 23.7, 22.3, 23.8, 23.1, 25.8, 25.5, 32.8, 33…
## $ af <dbl> 22.9, 23.6, 22.2, 24.4, 24.6, 19.1, 21.7, 21.1, 18.9, 19.8, 19…
## $ amf <dbl> 30.0, 28.4, 26.9, 26.7, 26.9, 27.1, 26.5, 24.7, 25.4, 21.7, 19…
## $ silte <dbl> 1.2, 1.2, 1.2, 0.6, 2.1, 2.2, 0.7, 0.2, 2.5, 0.2, 2.5, 2.6, 0.…
## $ argila <dbl> 21.5, 20.4, 21.4, 20.5, 20.7, 23.5, 23.1, 22.7, 22.0, 18.5, 20…
## $ s_a <dbl> 0.05, 0.05, 0.05, 0.02, 0.10, 0.09, 0.03, 0.01, 0.11, 0.01, 0.…
## $ af_ag <dbl> 6.16, 5.53, 4.44, 6.42, 7.94, 5.03, 6.03, 4.49, 4.20, 3.36, 3.…
## $ p <dbl> 42, 22, 41, 27, 11, 12, 11, 16, 38, 25, 25, 6, 6, 7, 5, 4, 4, …
## $ p_h <dbl> 4.2, 3.8, 4.8, 4.0, 4.4, 4.0, 4.8, 5.4, 4.4, 5.2, 4.5, 5.1, 4.…
## $ k <dbl> 0.27, 0.11, 0.34, 0.13, 0.11, 0.14, 0.23, 0.28, 0.19, 0.14, 0.…
## $ ca <dbl> 1.4, 0.4, 2.4, 0.7, 1.4, 0.6, 1.6, 3.3, 1.6, 2.9, 1.3, 1.6, 0.…
## $ mg <dbl> 0.3, 0.1, 0.4, 0.1, 0.3, 0.1, 0.7, 1.3, 0.5, 1.7, 0.6, 0.8, 0.…
## $ h_al <dbl> 5.2, 5.8, 4.2, 5.2, 4.2, 5.2, 3.4, 2.5, 5.2, 3.1, 4.2, 2.5, 4.…
## $ sb <dbl> 1.97, 0.61, 3.14, 0.93, 1.81, 0.84, 2.53, 4.88, 2.29, 4.74, 1.…
## $ t <dbl> 7.17, 6.41, 7.34, 6.13, 6.01, 6.04, 5.93, 7.38, 7.49, 7.84, 6.…
## $ v <dbl> 27, 10, 43, 15, 30, 14, 43, 66, 31, 60, 32, 50, 22, 35, 36, 34…dados |>
ggplot(
aes(
x=x,
y=argila,
color=sup
))+
geom_point()
#Aplicando um moelo de delineamento inteiramente casualizado
y<-dados |> pull(argila)
trat<-dados |> pull(sup) |>as_factor()
trat## [1] I I I I I I I I I I I I I I I I I II
## [19] II II II II II II II II II II II II II II II II II II
## [37] II II II II II II II II II II II II II II II II II II
## [55] II II II II II II II II II II II II II II II II II II
## [73] II II II II II II II III III III III III III III III III III III
## [91] III III III III III III III III III III III III III III III III
## Levels: I II III##Criar modelo para ánalise
mod<-aov(y~trat)
mod## Call:
## aov(formula = y ~ trat)
##
## Terms:
## trat Residuals
## Sum of Squares 1255.866 2727.294
## Deg. of Freedom 2 103
##
## Residual standard error: 5.145734
## Estimated effects may be unbalanced##Estrutura do objeto mod
str(mod)## List of 13
## $ coefficients : Named num [1:3] 21.1 -6.34 -10.96
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "tratII" "tratIII"
## $ residuals : Named num [1:106] 0.4 -0.7 0.3 -0.6 -0.4 ...
## ..- attr(*, "names")= chr [1:106] "1" "2" "3" "4" ...
## $ effects : Named num [1:106] -150.336 1.948 -35.385 -0.611 -0.411 ...
## ..- attr(*, "names")= chr [1:106] "(Intercept)" "tratII" "tratIII" "" ...
## $ rank : int 3
## $ fitted.values: Named num [1:106] 21.1 21.1 21.1 21.1 21.1 ...
## ..- attr(*, "names")= chr [1:106] "1" "2" "3" "4" ...
## $ assign : int [1:3] 0 1 1
## $ qr :List of 5
## ..$ qr : num [1:106, 1:3] -10.2956 0.0971 0.0971 0.0971 0.0971 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:106] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "tratII" "tratIII"
## .. ..- attr(*, "assign")= int [1:3] 0 1 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ trat: chr "contr.treatment"
## ..$ qraux: num [1:3] 1.1 1.11 1.16
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-07
## ..$ rank : int 3
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 103
## $ contrasts :List of 1
## ..$ trat: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ trat: chr [1:3] "I" "II" "III"
## $ call : language aov(formula = y ~ trat)
## $ terms :Classes 'terms', 'formula' language y ~ trat
## .. ..- attr(*, "variables")= language list(y, trat)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "trat"
## .. .. .. ..$ : chr "trat"
## .. ..- attr(*, "term.labels")= chr "trat"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, trat)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "trat"
## $ model :'data.frame': 106 obs. of 2 variables:
## ..$ y : num [1:106] 21.5 20.4 21.4 20.5 20.7 23.5 23.1 22.7 22 18.5 ...
## ..$ trat: Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ trat
## .. .. ..- attr(*, "variables")= language list(y, trat)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "trat"
## .. .. .. .. ..$ : chr "trat"
## .. .. ..- attr(*, "term.labels")= chr "trat"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, trat)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "trat"
## - attr(*, "class")= chr [1:2] "aov" "lm"##Extraindo os residuos do modelo
rs<-rstudent(mod)
rs## 1 2 3 4 5 6
## 0.079739246 -0.139552651 0.059803619 -0.119613529 -0.079739246 0.478958251
## 7 8 9 10 11 12
## 0.398994804 0.319106208 0.179436025 -0.518972767 -0.119613529 -0.379016227
## 13 14 15 16 17 18
## -0.199380742 0.458959546 -0.019934229 -0.099675805 -0.299143988 1.425307464
## 19 20 21 22 23 24
## 1.225272340 1.125799089 0.730837666 2.140207705 2.370928509 2.370928509
## 25 26 27 28 29 30
## 1.305103402 1.912650126 1.165548006 1.305103402 1.425307464 0.750492052
## 31 32 33 34 35 36
## 1.205350446 1.425307464 1.125799089 0.066039221 0.848895796 0.534710709
## 37 38 39 40 41 42
## 1.185442416 0.554292664 1.205350446 0.691553436 0.730837666 0.750492052
## 43 44 45 46 47 48
## 1.974419977 0.789826504 0.202556676 -0.148443870 -0.813952274 -0.499983185
## 49 50 51 52 53 54
## 0.007547180 -0.031446728 0.495564613 -0.519554627 -0.715627048 -1.110426371
## 55 56 57 58 59 60
## -1.389703543 -1.249711993 -1.429842586 -1.269665614 -0.774596123 -1.289633882
## 61 62 63 64 65 66
## -1.651851288 -0.774596123 -1.530485712 -1.570866517 -1.733147594 -0.912507996
## 67 68 69 70 71 72
## -1.189936578 -1.031119648 -0.187453873 -1.170039114 -0.813952274 -0.480417361
## 73 74 75 76 77 78
## -0.480417361 -1.409764874 -1.530485712 -1.209847718 -1.490176533 -1.189936578
## 79 80 81 82 83 84
## -1.329615342 -0.403227845 -0.265046270 -0.978988707 -0.324244221 0.089781704
## 85 86 87 88 89 90
## -0.919081451 -0.978988707 -1.159352014 -0.978988707 -0.939039543 -0.879197047
## 91 92 93 94 95 96
## -0.008758848 -0.186157061 -0.205875607 1.895996678 0.464709923 0.267238263
## 97 98 99 100 101 102
## 0.267238263 0.981209424 1.444431251 0.188347901 0.444942530 0.781869685
## 103 104 105 106
## -0.700194582 0.484482647 0.801759078 0.841565924##Extraindo os preditos pelo modelo
yp<-predict(mod)
yp## 1 2 3 4 5 6 7 8
## 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000
## 9 10 11 12 13 14 15 16
## 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000 21.10000
## 17 18 19 20 21 22 23 24
## 21.10000 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 25 26 27 28 29 30 31 32
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 33 34 35 36 37 38 39 40
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 41 42 43 44 45 46 47 48
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 49 50 51 52 53 54 55 56
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 57 58 59 60 61 62 63 64
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 65 66 67 68 69 70 71 72
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129
## 73 74 75 76 77 78 79 80
## 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 14.76129 10.14444
## 81 82 83 84 85 86 87 88
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444
## 89 90 91 92 93 94 95 96
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444
## 97 98 99 100 101 102 103 104
## 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444 10.14444
## 105 106
## 10.14444 10.14444Começar o diagnostico construindo um arquivo cm sup,y,yp,rs
diagnostico<-tibble(
trat,
y,
yp,
rs
)
diagnostico## # A tibble: 106 × 4
## trat y yp rs
## <fct> <dbl> <dbl> <dbl>
## 1 I 21.5 21.1 0.0797
## 2 I 20.4 21.1 -0.140
## 3 I 21.4 21.1 0.0598
## 4 I 20.5 21.1 -0.120
## 5 I 20.7 21.1 -0.0797
## 6 I 23.5 21.1 0.479
## 7 I 23.1 21.1 0.399
## 8 I 22.7 21.1 0.319
## 9 I 22 21.1 0.179
## 10 I 18.5 21.1 -0.519
## # ℹ 96 more rows#Normalidade dos residuos
diagnostico |>
ggplot(
aes(
x=rs
)
) +
geom_histogram(
bins=10,
color="black",
fill="gray30"
)+
theme_bw()
#QQplot
diagnostico |>
ggplot(
aes(
sample=rs
)
)+
stat_qq()+
stat_qq_line(
color="blue"
)
diagnostico |>
ggplot(aes(sample = rs)) +
stat_qq() +
stat_qq_line(color = "blue") +
facet_wrap(~ trat) +
theme_bw()
library(plotly)##
## Anexando pacote: 'plotly'## O seguinte objeto é mascarado por 'package:ggplot2':
##
## last_plot## O seguinte objeto é mascarado por 'package:stats':
##
## filter## O seguinte objeto é mascarado por 'package:graphics':
##
## layout# Calcular quantis teóricos e empíricos
qq_data <- qqnorm(diagnostico$rs, plot.it = FALSE)
# Suponha que 'trat' seja um fator indicando grupos
# Converta em numérico apenas para ilustrar no 3D
group_numeric <- as.numeric(diagnostico$trat)
plot_ly(
x = ~qq_data$x,
y = ~qq_data$y,
z = ~group_numeric,
type = 'scatter3d',
mode = 'markers'
) %>%
layout(
scene = list(
xaxis = list(title = "Quantis Teóricos"),
yaxis = list(title = "Quantis Empíricos"),
zaxis = list(title = "Grupo")
)
)# install.packages("nortest")
#Aplicando os testes de normalidade
shapiro.test(rs)#nao rejeitamos H0##
## Shapiro-Wilk normality test
##
## data: rs
## W = 0.97323, p-value = 0.03052nortest::lillie.test(rs)#nao rejeitamos H0##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: rs
## D = 0.072002, p-value = 0.1941nortest::cvm.test(rs)#nao rejeitamos H0##
## Cramer-von Mises normality test
##
## data: rs
## W = 0.077302, p-value = 0.2222# Conclusão
## Os resíduos seguem uma distribuição NORMAL
##Gráfico dos 5 números ou Boxplot
diagnostico |>
ggplot(aes(
x=trat,
y=y,
fill=trat,
)
)+
geom_boxplot()
diagnostico |>
ggplot(aes(x = trat, y = y, fill = trat)) +
# Gráfico de violino com ajuste na largura
geom_violin(trim = FALSE, alpha = 0.7, scale = "width") +
# Adicionar pontos individuais com jitter para visualizar a distribuição
geom_jitter(width = 0.1, alpha = 0.5, size = 1.5, color = "black") +
# Adicionar a mediana de forma destacada
stat_summary(fun = median, geom = "crossbar",
width = 0.3, color = "white",
fatten = 2, linetype = "solid") +
# Ajuste de escala de cor (ex: viridis - requer library(viridis) ou scale_fill_viridis_c)
scale_fill_brewer(palette = "Pastel1") +
# Título e rótulos
labs(title = "Distribuição por Tratamento",
x = "Tratamento",
y = "Resposta",
fill = "Tratamento") +
# Tema minimalista
theme_minimal(base_size = 14) +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5, face = "bold")
)
#Teste de homocedasticidade
#Teste de Levene
lawstat::levene.test(y,trat)##
## Modified robust Brown-Forsythe Levene-type test based on the absolute
## deviations from the median
##
## data: y
## Test Statistic = 22.666, p-value = 6.96e-09#Bartlett
bartlett.test(y,trat)##
## Bartlett test of homogeneity of variances
##
## data: y and trat
## Bartlett's K-squared = 31.396, df = 2, p-value = 1.522e-07##estudar a regularidade
##calcular media e variancia dos valores de y em funçao dos trat
diagnostico |>
group_by(trat) |>
summarise(
log_media=log(mean(y)),
log_variancia=log(var(y))
) |>
ggplot(
aes(
x=log_media,
y=log_variancia
)
)+
geom_point()+
geom_smooth(
method = "lm",
se=FALSE
)## `geom_smooth()` using formula = 'y ~ x'
#Análise de regressão
modelo_linear<-lm(log_variancia~log_media,
data=diagnostico |>
group_by(trat) |>
summarise(
log_media=log(mean(y)),
log_variancia=log(var(y))
)
)#transformação Box and Cox